Intermediate Mass Ratio Black Hole Binaries: 
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We study black-hole binaries in the intermediate-mass-ratio regime 0.01 < q < 0.1 with a new 
technique that makes use of nonlinear numerical trajectories and efficient perturbative evolutions 
to compute waveforms at large radii for the leading and nonleading modes. As a proof-of-concept, 
we compute waveforms for q — 1/10. We discuss applications of these techniques for LIGO/ VIRGO 
data analysis and the possibility that our technique can be extended to produce accurate waveform 
templates from a modest number of fully-nonlinear numerical simulations. 
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Introduction: The dramatic breakthroughs in the numer- 
ical techniques to evolve black-hole binaries (BHB) [l|- 
0] transformed the field and made it possible to begin 
constructing waveform templates for use in LIGO and 
VIRGO gravitational wave data analysis and detection 
from highly-accurate, fully-nonlinear simulations. This 
promises to greatly aid in the detection and analysis of 
gravitational waves from astrophysical sources 0] ■ 

While LISA's sensitivity band is well suited to observe 
supermassivc BHBs, ground based observatories are more 
sensitive to BHBs with masses up to a few hundred so- 
lar masses. Both of these types of BHBs are most likely 
to have mass ratios in the range 1:10 - 1:100 This 
regime is hard to model with fully-nonlinear numerical 
simulations, which thus far have focused on the compa- 
rable mass regime. On the other hand, BHBs in the 
very small mass ratio regime can, in principle, be mod- 
eled accurately with perturbation theory, and significant 
theoretical effort was dedicated to the consistent com- 
putation of the self-force on a point-like source orbit- 
ing around Schwarzschild and Kerr black holes (see Q 
for a review of new developments in the field). These 
perturbative techniques allow one to formally compute 
geodesic deviations to second order. Although secular 
effects still need to be quantified, these techniques are 
expected to provide a good approximation for mass ra- 
tios q = mi/m 2 < 1/100. 

Fully-nonlinear BHB simulations with small mass- 
ratios are challenging and the smallest mass ratios pub- 
lished so far are q = 1/10 in the nonspinning case 0, H[ 
and q = 1/8 for highly-spinning BHBs [§]. Most fully- 
nonlinear numerical simulations published have been in 
the comparable mass regime, and consequently, most of 
the comparisons between NR and post-Newtonian (PN) 
waveforms have been in this regime (for comparisons of 



NR and PN waveforms from generic binaries, see 
and references therein). Similarly, phenomenological 
and EOBNR [l3| modeling of waveforms are based on 
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comparable-mass BHB simulations. The intermediate 
region of 1/100 < q < 1/10 provides a unique oppor- 
tunity to compare fully-nonlinear numerical and pertur- 
bative techniques. In this letter we introduce a novel ap- 
proach to deal with this intcrmcdiatc-mass-ratio regime. 
We perform a 'proof of principle' comparison, to deter- 
mine an upper bound on the error in our technique, be- 
tween fully-nonlinear numerical simulations and pertur- 
bative ones in the borderline case of a nonspinning BHB 
system with mass ratio q = 1/10. 

The key element introduced in our approach is the 
use of the numerically generated trajectories (in the 
moving puncture approach (HQ) °f a BHB during the 
late inspiral regime in a perturbative calculation of the 
waveform emitted by the binary. 

Full Numerical Simulations: In order to approxi- 
mate the late orbital parameters of a quasicircular 
inspiral we computed the momentum and puncture- 
position parameters using a resummed (particle limit) 
3.5 PN evolution of a q = 1/10, non-spinning binary 
from r = 50M to r = 7.25M. To compute the numerical 
initial data, we use the puncture approach fbj ] along 
with the TwoPunctures [HI thorn. We chose the 
puncture mass parameters such that the system had a 
total ADM mass M = 1 and mass ratio q = 1/10. 

We evolved this BHB data-set using the LazEv [16| 
implementation of the moving puncture approach 
Our code used the Cactus toolkit [l?} and the Carpet [l8| 
mesh refinement driver to provide a 'moving boxes' style 
mesh refinement. We used a modified gauge condition 
suggested in [l||. Our simulation used 11 levels of re- 
finement (around the smaller components), with central 
resolutions as high as Mj 368. 64, and 9 levels of refine- 
ment around the larger component. The outer bound- 
aries were located at 400M and the resolution in the 
boundary zone was h = 2.7777M for our finest resolu- 
tion run. The BHB performs four orbits prior to merger, 



2 



TABLE I: Initial puncture data, and final remnant parameters 



Xl 


6.60438 


X2 


-0.671518 


m l 


8.43895 x 10" 2 


ml 


0.907039 


Px 


-3.2671 x 10" 4 


Py 


4.04057 x 10" 2 


M h 


0.99600 


ah 


0.26081 


£ rad * 1000 


4.578 ± 0.684 


JL* * ioo 


3.420 ± 0.113 


Vkick[kms _1 ] 


67 ±35 


tcAH 


« 380M 



y 4 (h)-y/ 4 (h/l.2) 

(\|i„(h/1 .2)-i|f 4 (h/1 .2 ))*2.0736 

(V 4 (h/1 .2)-¥ 4 (W1 .2 Z ))*4.29982 



=. 0.00015 




FIG. 1: Convergence of the real part of the (£ — 2, m = 2) 
mode of 1/14 as computed by full numerical evolutions and 
extracted at Robs = 100M. Eighth-order convergence im- 
plies i>i(h) - tp 4 {h/1.2) = 4.29982(i/>4(/i/1.2) - ip 4 {h/1.2 2 )), 
while fourth-order convergence implies ip 4 (h) — ip 4 {h/1.2) = 
2.0736(t/) 4 (V1- 2 ) - ipi{h/l.2 2 )). Initially, the error in ip A is 
very small and dominated by grid noise, at t ~ 250M the 
fourth-order convergence of the ^4 algorithm is apparent, with 
the convergence quickly changing to eighth-order as trunca- 
tion errors in the metric and extrinsic curvature begin to dom- 
inate the errors in ip4. 



which occurs roughly 400M after the start of the simu- 
lation (see Table |T|). 

We found that the waveform exhibits eighth-order 
convergence during the late-inspiral (when the er- 
rors in the waveform arc dominated by phase errors 
associated with the orbital motion). The waveform 
exhibits noise which tends to mask the convergence 
at early times. As the simulation progresses, we first 
see fourth-order convergence of the Weyl scalar t/> 4 (we 
used a fourth-order algorithm to compute ipi) an d later 
eighth-order convergence, as the errors in the waveform 
due to truncation errors in the metric and extrinsic 
curvature become important. In Fig [T] we show the con- 
vergence of the waveform during the late-inspiral. Note 
the very good agreement of the two finest resolution runs. 

Perturbative evolutions: For the case of a nonspin- 
ning BHB in the small mass ratio regime, we can use 
a Schwarzschild black-hole (rather than Kerr) as the 



background space-time for our perturbation expansion. 
This allows us to use the simpler metric perturbation 
equations and the Regge- Wheeler- Zerilli (RWZ) formal- 
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where V, 



(even) / (odd) 



denotes the Zerilli and Regge- 



Wheeler potentials and 5^ on )/( odd ) j s ^he source term 
derived from the particle's trajectory, in the orbital plane 
[r p (i), $ p (i)]. The details are discussed in 22,23)]. and 
we use the definitions in [24| ■ The consistent use of higher 
order tracks in this equations has been discussed in [25[ 
(See Eq. (4.6) therein.) 

The coordinates used in numerical simulations are cho- 
sen to produce stable evolutions and correspond, ini- 
tially, to isotropic coordinates. Perturbative calculations, 
on the other hand, regularly make use of the standard 
Schwarzschild coordinates. The easiest way to relate 
the two is to translate the numerical tracks into the 
Schwarzschild coordinates used in Eq. (JTJ). This can be 
achieved by considering the late-time numerical coordi- 
nates that correspond to radial isotropic 'trumpet' sta- 
tionary 1 + log slices of the Schwarzschild spacetime (2(| . 
We obtain the explicit coordinate transformations follow- 
ing the procedure detailed in Ref. [27[. We perform the 
numerical integration of the perturbative equation ([1]) 
with Dirac's delta sources using the algorithm described 

in [mim. 

An interesting aspect of the source term in the pertur- 
bative evolution equations ([T]) is that 5^ en )/( odd ) _s, q 



as (1 — 2M/r p ) when the particle approaches the back- 
ground horizon (located at r = 2M). This leads to essen- 
tially free evolution of the Cauchy data once the particle 
'sits' on the horizon. In this sense perturbation theory 
naturally transitions to a 'close limit' evolution. This 
property means that once the black holes are close to 
each other, and form a common horizon, we no longer 
need the information of the full numerical track, since 
perturbation equations will describe the radiation essen- 
tially in terms of the universal quasi-normal ringing. In 
practice, we monitor the numerical tracks and switch off 
the source-terms when r p (t) ~ 2M. 

In Fig. [2l we plot the leading order (i = 2, m = 2) 
mode of the strain h calculated both using the numer- 
ical waveform and, perturbatively, from the numerical 
trajectories, while in Fig. |3, we compare the nonlcading 
(£ = 2,m = 1) and (I = 3, m = 3) modes. The agree- 
ment between waveforms is summarized in TablellTlwhere 
we computed the overlap between waveforms as defined 
in Ref. [13]. In addition we evaluate the effects of ex- 
traction of the waveform at finite radii by extracting the 
perturbative waveform at i? b s = 100M and 1000M and 
then (after a time shift) computing the overlap index (see 
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FIG. 2: The Real part of the (£ = 2, m = 2) mode of the 
strain /i as computed by full numerical evolutions and ex- 
tracted at R = 100M, and the corresponding perturbative 
evolutions for the same numerical track extracted at two radii 
Robs ~ 100M, lOOOAf. The matching of the numerical and 
perturbative waveforms (i? bs = 100M and lOOOM) for the 
whole range of evolution for 130.068 < t/M < 593.564 are 
0.987656 and 0.974751, respectively; while for the two pertur- 
bative waveforms at different extraction radii the matching is 
0.938873. 



FIG. 3: Upper: The Real part of the (£ = 2, m = 1) mode 
of the strain h as computed by a numerical evolution and 
extracted at a radius i? b s = 100M, and the corresponding 
perturbative (odd parity) evolution for the same numerical 
track extracted at the same radius. Lower: The Real part of 
the [i = 3, m = 3) mode of h as computed by the numerical 
and perturbative evolutions. Similar results are obtained for 
the imaginary part of the modes. 
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TABLE II: The overlap (matching) of the real and imaginary 
parts of the modes of the strain h. The integration time 
is from t/M = 130.068 to 593.564 and the definition of the 
matching is given in Eqs. (26) and (27) of [l(J. 
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Match 


0.987656 




0.975017 




0.954138 
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1) 


= 3,771 = 


3) 
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0.987622 




0.981982 




0.953767 





Table [TTJ) . The equivalent extraction using standard nu- 
merical methods requires large computational resources, 
and such extractions have only recently been achieved 
using highly- efficient techniques such as multi-patch [29[ 
and pseudo-spectral methods. 

As we mentioned, the perturbative evolution only 
makes use of numerical information (tracks) in the 
final inspiral and merger phase. After merger we have 
free Cauchy evolution of the perturbations. One can 
extend this technique to the early inspiral phase, as 
well. During the early-inspiral, the binary's evolution is 
well described by PN theory and these trajectories can 
be used to determine the perturbative waveforms f3lj ) . 
We can thus, stitch together a 3.5PN trajectory for the 
inspiral phase to a full numerical evolution in the moving 
puncture approach. In fact, this is how we determined 
the initial orbital parameters for the numerical solution 
of General Relativity constraints to prepare full numeri- 
cal evolution. The resulting waveform is shown in Fig. 21 
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FIG. 4: The strain computed by a perturbation theory cal- 
culation of using the tracks provided by 3.5PN theory in the 
early-inspiral, the numerical tracks during the late-inspiral, 
and an effective close-limit calculation post-merger. Note that 
we obtain a waveform of duration ss 1600M using only the in- 
formation from the first 373M of a fully-nonlinear numerical 
evolution. 



which also shows the waveform produced solely with PN 
methods. Note that this procedure allows us to compute 
the complete waveform making use of numerical data 
from the first t « 373A/ (of an over 600M simulation) of 
the numerical evolution. The procedure resembles the 
'Lazarus approach' [32j in that it makes efficient use of 
numerical simulations by restricting them to the stage 
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of the binary evolution where they are needed. 

Conclusions and Discussion: Given the current dif- 
ficulties in efficiently solving the intermediate mass 
ratio regime by purely nonlinear numerical methods, we 
propose to use perturbation theory to propagate gravita- 
tional waves on a black hole background using trajectory 
information from a combination of nonlinear numerical 
and semi-analytic methods, such as the post-Newtonian 
approximation. In this letter we provided a 'proof of 
principle' comparison between our hybrid technique and 
fully nonlinear results in the most challenging case for 
perturbation theory, i.e. q = 1/10. 

The benefits of this approach are apparent: i) No 
fitting parameters have been used in the generation of 
the waveforms, ii) Since perturbation equations natu- 
rally transition to close limit ones, i.e. source free wave 
equations, numerical tracks are not needed after merger, 
thus reducing the cost of the nonlinear simulations (in 
our example it represented a 40% reduction in running 
time.) iii) The perturbative evolution also provides a 
natural way to extract waveforms at very large radii in a 
well known and defined gauge, in a region accurately de- 
scribed by linear perturbations. Thus further cutting the 
size of the required numerical grid and hence the compu- 
tational cost of the nonlinear simulation, iv) The com- 
putation of even, odd, leading, and higher-order modes 
(£, m) with the same black holes trajectory. Note that 
the small amplitude of the non-leading terms impose ad- 
ditional requirements of accuracy at the extraction radii 
and at the external boundaries of the full numerical grid 
to avoid interfering reflections, v) Perturbation codes 
are highly efficient and run within seconds on standard 
CPUs, and are also amenable for running in GPUs or 
other accelerated hardware opening up the possibility of 
generating online analysis of LIGO/ VIRGO data. 

In order to improve the accuracy of the generated 
waveforms we can consider several extensions to our 
method. These include perturbations about a Kerr back- 
ground using the Teukolsky equations [33| . In addition, 
adding second-order perturbative corrections [34j should 
improve accuracy for larger mass ratios. Finally, given 
our encouraging first results, one can consider a phe- 
nomenological description of intermediate-mass-ratio in- 
spiral trajectories based on a series of detailed nonlincar 
numerical runs to model the tracks of BHBs as function 
of mass ratio and spins. The resulting phenomenologi- 
cal trajectories could then be used to efficiently generate 
accurate waveforms over a wide region of the parameter 
space. 
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